---
title: "JLM_CN_code"
author: "Jonathan Messerschmidt"
date: "2025-11-05"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
# Load libraries
library(Seurat)
library(tidyverse)
# library(biomaRt)
library(dplyr)
library(gridExtra)
library(patchwork)
library(ggplot2)
library(grid)
library(data.table)
library(RColorBrewer)
# library(monocle3)
library(Matrix)

library(readxl)
library(ggpubr)
library("R.utils")
# library("biomaRt")
# library(AnnotationDbi)
# library(msigdbr)
# library(fgsea)

library(svglite)

# 8 / 15 / 24
library(broom)
library(pheatmap)
library(SCpubr)

# 8/28/24
library(anticlust)
library(factoextra)
library(cluster)
source("clustering-functions.R")
# install.packages('devtools') #### If needed to install the presto package to enhance DEG determination
# devtools::install_github('immunogenomics/presto')

library(here)

```

# Load Seurat Objects
```{r}
seu.file <- "seu_obj_cc_020924.rds"
seu.obj.cc <- readRDS(seu.file)
```

# Small QC plots
```{r}
seu.obj.cc@meta.data$all <- NULL
seu.obj.cc@meta.data$all <- "all"

Idents(seu.obj.cc) <- "all"
p <- VlnPlot(seu.obj.cc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), alpha = 0)


Idents(seu.obj.cc) <- "treatment"

plotlist <- list(
  VlnPlot(seu.obj.cc, idents = "WT", features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), alpha = 0),
  VlnPlot(seu.obj.cc, idents = "OPNKO", features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), alpha = 0),
  p
  )


cowplot::plot_grid(plotlist = plotlist, nrow = 3)
```


# Obtain the cell IDs of astrocytes and microglia for subsetting.
8/15/24: Edited code to obtain T cells instead of Endothelial cells
10/7/25: Reviewer Comment --> dotplot vs. heatmap
```{r}
Idents(seu.obj.cc) <- "celltype"
DefaultAssay(seu.obj.cc) <- "RNA"

astrocyte_IDs.WT <- rownames(seu.obj.cc@meta.data[
  seu.obj.cc@meta.data$celltype == "Astrocytes" & seu.obj.cc@meta.data$treatment == "WT",
  ])

microglia_IDs.WT <- rownames(seu.obj.cc@meta.data[
  seu.obj.cc@meta.data$celltype == "Microglia" & seu.obj.cc@meta.data$treatment == "WT",
  ])

oligo_IDs.WT <- rownames(seu.obj.cc@meta.data[
  seu.obj.cc@meta.data$celltype == "Oligodendrocytes" & seu.obj.cc@meta.data$treatment == "WT",
  ])

T_IDs.WT <- rownames(seu.obj.cc@meta.data[
  seu.obj.cc@meta.data$celltype == "T Cells" & seu.obj.cc@meta.data$treatment == "WT",
  ])
```

# Subset the data based on cell IDs, normalize and obtain DEGs.
```{r}
## subsetting
astro <- subset(seu.obj.cc, cells = astrocyte_IDs.WT) 
micro <- subset(seu.obj.cc, cells = microglia_IDs.WT)
oligo <- subset(seu.obj.cc, cells = oligo_IDs.WT) 
tcells <- subset(seu.obj.cc, cells = T_IDs.WT) 

# Scale and Normalize the RNA data layer

obj.list <- list(astro = astro, micro = micro, oligo = oligo, tcells = tcells)

obj.list <- lapply(X = obj.list, FUN = function(obj){
  
  obj <- obj %>% 
    NormalizeData() %>% 
    FindVariableFeatures(nfeatures = 2000) %>% 
    ScaleData()
  
  return(obj)
})



# This portion reorders the time.point factor to make more sense of it
obj.list <- lapply(obj.list, function(obj){
  
  Idents(obj) <- "time.point"
  time.ordered.obj <- obj@meta.data$time.point
  time.ordered.obj <- factor(time.ordered.obj, levels = c("0d", "4d", "14d"))
  obj$time.point <- time.ordered.obj
  Idents(obj) <- "time.point"
  
  return(obj)
})
```

# Portion below for Figure of Microglial 0, 4, 14 DEGs --> 10/7/25 Redoing as DotPlot
```{r}
##########
# Checking to make sure that the two cluster-like metadata features are the same (TRUE)
identical(obj.list$micro$integrated_snn_res.0.7, obj.list$micro$seurat_clusters)

SCpubr::do_DimPlot(obj.list$micro, group.by = "seurat_clusters")
SCpubr::do_DimPlot(obj.list$micro, group.by = "time.point")

# 10/13/25 --> Adding Astrocyte dotPlots
SCpubr::do_DimPlot(obj.list$astro, group.by = "seurat_clusters")
SCpubr::do_DimPlot(obj.list$astro, group.by = "time.point")

# UMAP coordinates have not been recalculated --> proceed to silhouette-optimized resolution clustering

#
#
# Create a new object that can then be reclustered
# Change out obj.list$micro for the object of interest

############### BEGIN REPEAT FOR EACH OBJECT AS NEEDED FOR UMAP / CLUSTERING SOLUTION #################### 

obj.normalized <- obj.list$astro
#
#

obj.normalized <- obj.normalized %>% 
  RunPCA(features = VariableFeatures(.))

# Take a look at how the PCs split the data
ElbowPlot(obj.normalized)
# Top 15 PCs carried for further analysis

# Top 15 features per PC
VizDimLoadings(obj.normalized, dims = 1:5, nfeatures = 15, col = "blue")

# Plot the PCA
SCpubr::do_DimPlot(obj.normalized, reduction = "pca", dims = c(1,2))

# Run the UMAP
obj.normalized <- RunUMAP(obj.normalized, dims = 1:15, metric = "correlation")
SCpubr::do_DimPlot(obj.normalized, reduction = "umap")

obj.normalized <- FindNeighbors(obj.normalized, dims = 1:15)
obj.normalized <- FindClusters(obj.normalized, resolution = seq(0.1, 1, by = 0.05))

top15_PCs <- as.data.frame(obj.normalized$pca@cell.embeddings[,1:15])
euc_distance <- dist(top15_PCs, method = "euclidean")

Louvain_SilScores <- Louvain_SilhouetteScore_RNA(obj.normalized, euc_distance)

#### Plot the average Silhouette Score values 
ggplot(data.frame(x = seq(0.1, 1, by = 0.05), y = Louvain_SilScores), 
       aes(x = x, y = y)) + 
  geom_point(pch = 18, size = 3) + 
  labs(x = "Resolution Parameter", y = "Silhouette Scores", 
       title = "Microglia Cluster Louvain Average Silhouette Scores") + # Change cell type in plot
  theme(axis.text.x = element_text(angle = 90))

# At this point, you have to choose which resolution to use for cluster identification (maximal Silhouette Score).
# In this case (microglia), the minimal resolution (0.1) is chosen.

fviz_silhouette(silhouette(as.numeric(obj.normalized$RNA_snn_res.0.1), euc_distance),
                palette=c("salmon", "darkblue", "dodgerblue"))

#astrocytes
fviz_silhouette(silhouette(as.numeric(obj.normalized$RNA_snn_res.0.2), euc_distance),
                palette=brewer.pal(8, "Blues"))

# Assigning the identity of the clusters to the cells and visualizing solution.
# Micro == 0.1 resolution chosen
# Astro == 0.2 resolution chosen
Idents(obj.normalized) <- "RNA_snn_res.0.2"
SCpubr::do_DimPlot(obj.normalized, group.by = "RNA_snn_res.0.2") # Optimal clustering solution
SCpubr::do_DimPlot(obj.normalized, group.by = "time.point")

micro.normalized <- obj.normalized
astro.normalized <- obj.normalized

############### END OF REPEAT FOR CLUSTERING SOLUTION #####################
SCpubr::do_DimPlot(micro.normalized, group.by = "RNA_snn_res.0.1") # Optimal clustering solution
SCpubr::do_DimPlot(micro.normalized, group.by = "time.point")
SCpubr::do_DimPlot(astro.normalized, group.by = "RNA_snn_res.0.2") # Optimal clustering solution
SCpubr::do_DimPlot(astro.normalized, group.by = "time.point")


# Find markers differentially expressed between identities of choice.
Idents(micro.normalized) <- "RNA_snn_res.0.1"
micro.normalized.markers <- FindMarkers(micro.normalized, ident.1 = "1", ident.2 = c("0", "2"), logfc.threshold = 0.25,
                             only.pos = TRUE)

# 10/7/25 --> Top 20 markers
# Plotting Markers of interest
top.celltype <- micro.normalized.markers %>% rownames_to_column(var = "gene_name") %>%
  arrange(p_val_adj) %>% dplyr::slice(1:20) %>% pull(gene_name) %>% rev()


micro.normalized$time.point <- factor(micro.normalized$time.point, levels = c("14d", "4d", "0d"))

p.yl.rd.no_black <- SCpubr::do_DotPlot(micro.normalized, group.by = "time.point", flip = T,
                   features = top.celltype,
                   legend.position = "right", dot.scale = 8)+
  scale_fill_gradient2(low = "white", mid = "orange2", high = "red2", midpoint = 2)


micro.genes <- c("Cx3cr1", "Gpr34", "P2ry12", "Hexb", "Tmem119")
p.microGenes.yl.rd.no_black <- SCpubr::do_DotPlot(micro.normalized, group.by = "time.point", flip = T,
                   features = micro.genes,
                   legend.position = "right", dot.scale = 8)+
  scale_fill_gradient2(low = "white", mid = "orange2", high = "red2", midpoint = 2.5)


#### ADDITIONAL DOT PLOTS -- MG GENES , REACTIVE ASTROCYTE SIGNATURE ########################

astro.normalized$time.point <- factor(astro.normalized$time.point, levels = c("14d", "4d", "0d"))

features.astro <- c("Gbp2", "Cxcl10", "Stat3","Osmr",  "Vim", "Lcn2", "Gfap", "Cd44")

p.astro_sig.yl.rd.no_black <- SCpubr::do_DotPlot(astro.normalized, group.by = "time.point", flip = T,
                   features = features.astro,
                   legend.position = "right", dot.scale = 8)+
  scale_fill_gradient2(low = "white", mid = "orange2", high = "red2", midpoint = 0.5)
```


# Noted differences in Expression --> UMAP
```{r}
# Run above to generate micro.normalized and astro.normalized
Idents(micro.normalized) <- "time.point"
p1 <- SCpubr::do_DimPlot(micro.normalized, plot.title = "Microglia, All")

p2 <- SCpubr::do_FeaturePlot(micro.normalized,
                       features = "Tmem119")
p3 <- SCpubr::do_FeaturePlot(micro.normalized,
                       features = "P2ry12")
p4 <- SCpubr::do_FeaturePlot(micro.normalized,
                       features = "Gpr34")
p5 <- SCpubr::do_FeaturePlot(micro.normalized,
                       features = "Hexb")
p6 <- SCpubr::do_FeaturePlot(micro.normalized,
                       features = "Cx3cr1")

layout <- "
##BBCC
AADDEE
####FF
"
p <- p1 + p2 + p3 + p4 + p5 + p6 +
  plot_layout(design = layout)

```


```{r}
# Back to original code below...
# Here I find the markers differentially expressed between 0d and 14d, then filtered.
obj.list.markers <- lapply(obj.list, function(obj){
  
  obj.markers <- FindMarkers(obj, ident.1 = "14d", ident.2 = "0d", logfc.threshold = 0.1,
                             only.pos = TRUE)
  
  obj.markers.filt <- obj.markers %>% 
    filter(p_val_adj < 0.05 & abs(avg_log2FC) > 0.25)
  
  return(list(raw = obj.markers, filt = obj.markers.filt))
})

obj.list.markers.raw <- list(astro = obj.list.markers$astro$raw,
                             micro = obj.list.markers$micro$raw,
                             oligo = obj.list.markers$oligo$raw,
                             tcells = obj.list.markers$tcells$raw)
```

# Below here is where GSEA was performed --- Not in manuscript


__________

# IFN-gamma DEG heatmaps for MG and astrocytes

---
title: "IFNy-DEGs"
author: "Jon Messerschmidt"
date: "2024-04-03"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
# Load libraries
library(Seurat)
library(tidyverse)
library(biomaRt)
library(dplyr)
library(gridExtra)
library(patchwork)
library(ggplot2)
library(grid)
library(data.table)
library(RColorBrewer)
# library(monocle3)
library(Matrix)
```

# Load Seurat Objects
```{r}
seu.file <- "seu_obj_cc_020924.rds"
seu.obj.cc <- readRDS(seu.file)
```

# Obtain the cell IDs of astrocytes and microglia for subsetting, then read in the IFNy hallmark geneset.
```{r}
Idents(seu.obj.cc) <- "celltype"
astrocyte_IDs.WT <- rownames(seu.obj.cc@meta.data[
  seu.obj.cc@meta.data$celltype == "Astrocytes" & seu.obj.cc@meta.data$treatment == "WT",
  ])

microglia_IDs.WT <- rownames(seu.obj.cc@meta.data[
  seu.obj.cc@meta.data$celltype == "Microglia" & seu.obj.cc@meta.data$treatment == "WT",
  ])

ifny <- read.delim("HALLMARK_INTERFERON_GAMMA_RESPONSE.v2023.2.Mm.tsv", sep = "\t")
ifny.genes <- ifny %>% filter(STANDARD_NAME == "GENE_SYMBOLS") %>% 
  pull(HALLMARK_INTERFERON_GAMMA_RESPONSE) %>% 
  strsplit(., ",") %>% 
  unlist()

DefaultAssay(seu.obj.cc) <- "RNA"

### 4/9/24 Request for oligodendrocytes and endothelial cells
Idents(seu.obj.cc) <- "celltype"
oligo_IDs.WT <- rownames(seu.obj.cc@meta.data[
  seu.obj.cc@meta.data$celltype == "Oligodendrocytes" & seu.obj.cc@meta.data$treatment == "WT",
  ])

endo_IDs.WT <- rownames(seu.obj.cc@meta.data[
  seu.obj.cc@meta.data$celltype == "Endothelial Cells" & seu.obj.cc@meta.data$treatment == "WT",
  ])
```

# Subset the data based on astrocyte cell IDs and plot DEGs in IFNY hallmark pathway.
```{r}
## Astrocytes subsetting
astro <- subset(seu.obj.cc, cells = astrocyte_IDs.WT) 

# Scale and Normalize the RNA data layer
astro <- astro %>% 
  NormalizeData() %>% 
  FindVariableFeatures(nfeatures = 2000) %>% 
  ScaleData()

# This portion reorders the time.point factor to make more sense of it
Idents(astro) <- "time.point"
time.ordered.a <- astro@meta.data$time.point 
time.ordered.a <- factor(time.ordered.a, levels = c("0d", "4d", "14d"))
astro$time.point <- time.ordered.a
Idents(astro) <- "time.point"

# Here I find the markers differentially expressed between 0d and 14d. Note that this includes ALL DEGs filtered by the parameters below.
astro.markers <- FindMarkers(astro, ident.1 = "0d", ident.2 = "14d")
astro.markers <- astro.markers %>% 
  filter(p_val_adj < 0.05 & abs(avg_log2FC) > 0.25)

# Here I convert the genes to a vector and only keep those in the IFNY geneset
astro.genes <- as.vector(rownames(astro.markers))
features.a <- astro.genes[astro.genes %in% ifny.genes]

# Plotting the heatmap
DoHeatmap(astro, features = features.a, group.by = "time.point") + 
  scale_fill_gradientn(colors = c("blue", "white", "red")) + 
  labs(title = "Astrocytes: Significant IFNy Geneset DEGs (0d vs 14d)") +
  theme(plot.margin = margin(1.5,0.5,0,1, "cm"), plot.title = element_text(vjust = 5, hjust = 0.5))
```

# Same steps as above, but for microglia
```{r}
## Microglia Subsetting
micro <- subset(seu.obj.cc, cells = microglia_IDs.WT) 

# Renormalization
micro <- micro %>% 
  NormalizeData() %>% 
  FindVariableFeatures(nfeatures = 2000) %>% 
  ScaleData()

# Setting factor levels / order
Idents(micro) <- "time.point"
time.ordered.mg <- micro@meta.data$time.point 
time.ordered.mg <- factor(time.ordered.mg, levels = c("0d", "4d", "14d"))
micro$time.point <- time.ordered.mg
Idents(micro) <- "time.point"

# Finding DEGS
micro.markers <- FindMarkers(micro, ident.1 = "0d", ident.2 = "14d")
micro.markers <- micro.markers %>% 
  filter(p_val_adj < 0.05 & abs(avg_log2FC) > 0.25)

# Filtering DEGs for only IFNY hallmark genes
micro.genes <- as.vector(rownames(micro.markers))
features.mg <- micro.genes[micro.genes %in% ifny.genes]

# Plot
DoHeatmap(micro, features = features.mg, group.by = "time.point") + 
  scale_fill_gradientn(colors = c("blue", "white", "red")) + 
  labs(title = "Microglia: Significant IFNy Geneset DEGs (0d vs 14d)") +
  theme(plot.margin = margin(1.5,0.5,0,1, "cm"), plot.title = element_text(vjust = 5, hjust = 0.5))
```

